Coupling Lattice Boltzmann with Atomistic Dynamics for the multiscale simulation of 

nano-biological flows 
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We describe a recent multiscale approach based on the concurrent coupling of constrained molec- 
ular dynamics for long biomolecules with a mesoscopic lattice Boltzmann treatment of solvent hy- 
drodynamics. The multiscale approach is based on a simple scheme of exchange of space-time 
information between the atomistic and mesoscopic scales and is capable of describing self-consistent 
hydrodynamic effects on molecular motion at a computational cost which scales linearly with both 
solute size and solvent volume. For an application of our multiscale method, we consider the much 
studied problem of biopolymer translocation through nanopores: we find that the method repro- 
duces with remarkable accuracy the statistical scaling behavior of the translocation process and 
provides valuable insight into the cooperative aspects of biopolymer and hydrodynamic motion. 



I. INTRODUCTION 

Modeling biological systems in an efficient and reliable 
way is a delicate task, which calls for the development 
of innovative computational methods, often requiring so- 
phisticated upgrades and extensions of techniques origi- 
nally developed for physical and/or chemical stand-alone 
applications. Indeed, biological systems exhibit a de- 
gree of complexity and diversity straddling across many 
decades in space-time resolution, to the point that, for 
many years, biological systems served as a paradigm of 
the kind of complexity which can only be handled in 
qualitative or descriptive terms. Advances in computer 
technology, combined with constant progress and break- 
throughs in simulational methods, are closing the gap 
between quantitative models and actual biological be- 
havior. The main computational challenge raised by bi- 
ological systems remains the wide and disparate range of 
spatio-temporal scales involved in their dynamical evo- 
lution, with protein folding, morphogenesis, intra- and 
extra-cellular communication, being just a few examples. 

In response to this challenge, various strategies have 
been developed recently, which are in general referred to 
as multiscale modeling. These methods are based on com- 
posite computational schemes which rely upon multiple 
levels of description of a given biological system, most 
typically the atomistic and continuum levels. The multi- 
ple descriptions are then glued together, through suitable 
'hand-shaking' procedures, to produce the final compos- 
ite multiscale algorithm. To date, the mainstream multi- 
scale modeling is based on the coupling between atomistic 
and continuum models. This choice reflects the historical 
developments of statistical mechanics and computational 
physics. In essence, continuum methods reduce the in- 
formation to a small number of distributed properties 
(fields), whose space-time evolution is computed by solv- 
ing a corresponding set of partial differential equations, 
such as reaction-diffusion-advection equations. Atomistic 
models, on the other hand, rely on the well established 
approach of molecular dynamics, possibly including ex- 



tensions capable of dealing with a quantum description. 
It is perhaps interesting to notice that this two-stage 
continuum-to-atomistic representation overlooks a third, 
intermediate, level of description, that is the mesoscopic 
level, as typically represented by Boltzmann kinetic the- 
ory and its extensions. Kinetic theory lies between the 
continuum and atomistic descriptions, and it is thereby 
natural to expect that it should provide an appropri- 
ate framework for the development of robust multiscale 
methodologies. 

Until recently, this approach has been hindered by the 
fact that the central equation of kinetic theory, the Boltz- 
mann equation, was typically perceived computationally 
nearly as demanding as molecular dynamics, and yet of 
very limited use for dense fluids (water being the typical 
biological medium), due to the lack of many-body cor- 
relations. Recent developments in lattice kinetic theory 
[J0] are making this view obsolete. Over the last decade, 
such developments have provided solid evidence that suit- 
ably discretized forms of minimal kinetic equations, and 
most notably the Lattice Boltzmann equation, are giv- 
ing rise to very efficient algorithms capable of handling 
complex flowing systems across many scales of motion. 
The behavior of fluid flow is described through minimal 
forms of the Boltzmann equation, living on a discrete lat- 
tice. The lattice dynamics is designed in such a way as 
to reflect the basic conservation laws of continuum me- 
chanics, and also to host additional (mesoscopic) physics 
which is not easily accommodated by continuum mod- 
els. Remarkably, both tasks can be achieved within the 
same algorithm, which proves often computationally ad- 
vantageous over the continuum approach based on the 
Navier-Stokes equations. It is only very recently that LB 
advances have started to be incorporated within a new 
class of mesoscopic multiscale solvers [3| . In this work, we 
address precisely such mesoscopic multiscale solvers, with 
specific focus on the biologically important problem of 
biopolymer translocation through nanopores. Our proce- 
dure is based on the assumption that, in order to capture 
the essential aspects of the translocation process, it is not 
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necessary to resolve all underlying atomistic details. As 
a result, both solvent and solute degrees of freedom are 
treated through appropriate coarse-graining. Due to the 
intrinsic coarse-graining nature of our methodology, the 
direct mapping to experimental conditions has to proceed 
through the adjustment of appropriate parameters. 



II. MULTISCALE COUPLING METHODOLOGY 

We will discuss the implementation of how a meso- 
scopic fluid solver, the lattice Boltzmann method (LB), 
can be coupled concurrently to the atomistic scale em- 
ploying explicit atomistic dynamics which, for simplicity, 
will be named molecular dynamics (MD) in a broad sense. 
This procedure involves different levels of the statistical 
description of matter (continuum and atomistic) and is 
able to handle different scales through the spatial and 
temporal coupling between the constrained molecular dy- 
namics for the polymer evolution and the lattice Boltz- 
mann treatment of the explicit solvent dynamics. This 
multiscale framework is well suited to address a class of 
biologically related problems. 

The solvent dynamics does not require any form of sta- 
tistical ensemble averaging, as it is represented through 
a discrete set of pre-averaged probability distribution 
functions (the single-particle Boltzmann distributions), 
which are propagated along straight particle trajecto- 
ries. At variance with Brownian dynamics, the lattice 
Boltzmann approach handles the fluid-mediated solvent- 
solvent interactions through an explicit representation of 
local collisions between the solvent and solute molecules. 
By leveraging space-time locality, the corresponding al- 
gorithm scales linearly with the number of beads, as op- 
posed to the (super)quadratic dependence of Brownian 
dynamics. This dual field/particle nature greatly facili- 
tates the coupling between the mesoscopic and atomistic 
levels, both on conceptual and computational grounds. 
Full details on this scheme are reported in Ref. It 
should be noted that, LB and MD have been coupled be- 
fore for the investigation of single-polymer dynamics [5| , 
but here this coupling is extended to long molecules of 
biological interest. 

A word of caution is in order with the Stokes limit 
Re — > 0. In this limit, the scale separation between 
atomistic and hydrodynamic degrees of freedom becomes 
opaque, as the atomic mean free path becomes compa- 
rable with hydrodynamic scales. However, it has been 
shown that finite Reynolds corrections to hydrodynam- 
ics have very negligible effects on the solvent-mediated 
forces between suspended bodies, e.g. Oseen-level hy- 
drodynamics is satisfactorily recovered [1, H, 0, H] • 

We first turn to the atomistic part within our approach 
and consider a polymer consisting of N monomer units 
(also referred to as beads). Each bead of the polymer is 
advanced in time according to the following set of Molec- 



ular Dynamics equations: 

Ffot,i — * c,« ~t~ Fdrag.i F r i -\- F K i 



(1) 



where the index i runs over all beads. In this expression, 
F Cji is a conservative force describing bead-bead interac- 
tions, represented here by a Lennard-Jones potential: 



V LJ {r)=Ae 



(2) 



This potential is truncated at a distance of and 
augmented by an angular harmonic term to account for 
distortions of the angle between consecutive bonds: 
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with <j) the relative angle between two consecutive bonds, 
and K,p a constant. Torsional motions are not included 
in the present model, but can easily be incorporated if 
needed. The second term in Eq.(TT]), Fd ra g,i, represents 
the dissipative drag force due to polymer-fluid coupling 
given by 
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with Vi, Ui the bead and fluid velocity evaluated at the 
bead position r*j of bead i with a mass m; 7 is the friction 
coefficient. In addition to mechanical drag, the polymer 
feels the effects of stochastic fluctuations of the fluid en- 
vironment, which is related to the third term in Eq.(TT]), 
F r ^ , an uncorrelated random force with zero mean acting 
on bead i. The term F Kt i in Eq.([T]) is the reaction force 
resulting from N — 1 holonomic constraints for molecules 
modelled with rigid covalent bonds. The usage of con- 
straints instead of flexible bond lengths makes it possible 
to eliminate unimportant high-frequency intra-molecular 
motion which might lead to numerical instabilities. In or- 
der to avoid spurious dissipation, the bead velocities are 
required to be strictly orthogonal to the relative displace- 
ments. The constraints on both positions and velocities 
are enforced over positions and momenta sep arately via 
the SHAKE and RATTLE algorithms [1, [Of . 

The LB equation is a minimal form of the Boltzmann 
kinetic equation in which all details of molecular motion 
are removed except those that are strictly needed to re- 
cover hydrodynamic behavior at the macroscopic scale 
(mass-momentum and energy conservation). The result 
is an elegant equation for the discrete distribution func- 
tion f p (x,t) describing the probability to find a LB par- 
ticle at lattice site x at time t with a discrete speed c p . 
Specifically, in this work we are dealing with nanoscale 
flows and will consider the fluctuating Lattice Boltzmann 
equation which takes the following form [Til ]: 




MD time 



— — LB time 



(c) 

FIG. 1: Transfer of spatial information (a) from grid to particle, and (b) from particle to grid. Grey spheres denote beads, while 
in white are the lattice sites. In (c) the information exchange (LB-MD couling) is sketched through the vertical highlighted 
regions between the two different scales of our multiscale approach. The MD marches in time on a finer time-scale (by a factor 
of 5 in the sketch and our simulations) than the LB solver. 



f p (x + c p At, t + At) = f P (x, t) - wAt(/ p - f; q ){x , t) + FpAt + S p At 

I 



(5) 



The particles can only move along the links of a regular 
lattice denned by the discrete speeds, so that the syn- 
chronous displacements Ax p = c p At never take the fluid 
particles away from the lattice. For the present study, 
the standard three-dimensional 19-speed lattice is used 
(see Fig.l in Ref. 0). The right hand side of Eq. [5] rep- 
resents the effect of intermolecular solvent-solvent col- 
lisions, through a relaxation toward local equilibrium, 
fp q , typically a second order (low-Mach) expansion in 
the fluid velocity of a local Maxwellian with speed u: 



f; q = w pP {\ 



+ -^-a[uu : {CpC p - c 2 a *f)]} (6) 
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where c s is the sound speed of the solvent, w p is a set 
of weights normalized to unity, I is the unit tensor in 
configuration space, and p is the local density. The relax- 
ation frequency u> controls the fluid kinematic viscosity 
v, through the relation v = c 2 s (1/uj — At/2) || and At the 
LB time-step. Knowledge of the discrete distributions f p 
allows the calculation of the local density p, flow speed 
pu and momentum- flux tensor P , by a direct summation 
upon all discrete distributions: 



P(x,t) = ^2f p (x,t) 
p 

pu(x,t) = ^2fp(x,t)c p 

p 

P{x,t) = ^2fp(x,t)c p Cp 



(7) 
(8) 
(9) 



The diagonal component of the momentum-flux tensor 
gives the fluid pressure, while the off-diagonal terms give 



the shear-stress. Both quantities are available locally and 
at any point in the simulation. Thermal fluctuations are 
included through the source term F p in Eq. l[5]). which 
is consistent with the fluctuation-dissipation theorem at 
all scales. In the same equation, the polymer-fluid back 
reaction is described through the source term S p , which 
represents the momentum input per unit time due to the 
reaction of the polymer on the fluid populations. This 
back reaction is given by the following expression: 



Sp (x,t) 
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s i€D(x) 



drag,i 



(10) 



where D{x) denotes the mesh cell to which the t h bead 
belongs. All quantities in this equations have to reside 
on the lattice nodes, thereby the frictional and random 
forces need to be extrapolated from the particle to the 
grid location. 

In the LB solver, free-streaming proceeds along 
straight trajectories which secures exact conservation of 
mass and momentum of the numerical scheme, but also 
greatly facilitates the imposition of geometrically com- 
plex boundary conditions. There is no need to solve the 
computationally expensive Poisson equations, since the 
pressure field is locally available. All interactions are 
local, rendering the LB scheme ideal for parallel comput- 
ing. More advanced Lattice Boltzmann models [l2j also 
have been developed and could equally well be suited for 
coupling to atomic scale dynamics. 

The Molecular Dynamics solver is marched in 
time with a stochastic integrator (due to extra non- 
conservative and random terms) [131 ], proceeding at a 
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fraction 1/M of the LB time-step At: At M D = At/M. 
The time-step ratio M > 1 controls the scale separation 
between the solvent and solute timescales and should be 
chosen as small as possible, consistent with the require- 
ment of providing a realistic description of the polymer 
dynamics. The MD cycle is repeated M times, with 
the hydrodynamic field frozen at each LB time-stamp 
t n = nAt, at which the transfer of spatial information 
from grid to particle locations (and conversely) is per- 
formed. For this transfer, a simple nearest grid point 
interpolation scheme is used (see FigfTJa), (b)), on ac- 
count of its simplicity. At a time step t n = nAt, the 
pseudo-algorithm describing single LB time-step, will be: 

1. Interpolation of the velocity: u(x) — > Ui 

2. For k = 1, M : Advance the molecular state from t 
to t + AtMD 

3. Extrapolation of the forces: Fi — > F(x) 

4. Advance the Boltzmann populations from t to t + 
At 

A sketch of this scheme is presented in Fig. [TJ In terms 
of computational efficiency, At m d is largely independent 
of the number of beads N because (i) the LB-MD cou- 
pling is local, (ii) the forces are short ranged and (iii) the 
SHAKE/RATTLE algorithms are empirically known to 
scale linearly with the number of constraints. Up to this 
point we described a scheme that is general and applica- 
ble to any situation where a long polymer is moving in a 
solvent. This motion is of great interest for a fundamen- 
tal understanding of polymer dynamics in the presence of 
the solvent and crucial in relevant biophysical processes. 



III. BIOPOLYMER TRANSLOCATION 
THROUGH NANOPORES 

We next turn to an application of the multiscale 
scheme described in the previous section. We were mo- 
tivated by recent experimental studies which have fo- 
cused on the translocation of biopolymers such as RNA 
or DNA through nanometer sized pores. These explore in 
vitro the translocation process through micro-fabricated 
channels under the effects of an external electric field, or 
through protein channels across cellular membranes [lj] . 
In particular, recent experimental work has focussed on 
the possibility of fast DNA-sequencing through electronic 
means, that is, by reading the DNA sequence while it is 
moving through a nanopore under the effect of a localized 
electric field [lj] • This type of biophysical processes are 
important in phenomena like viral infection by phages, 
inter-bacterial DNA transduction or gene therapy [la ]. 
Some universal features of DNA translocation have also 
been analyzed theoretically, by means of suitably sim- 
plified statistical schemes Jig], and non-hydrodynamic 
coarse-grained or microscopic models fl7l fislj . However, 



these complex phenomena involve the competition be- 
tween many-body interactions at the atomic or molecu- 
lar scale, fluid-atom hydrodynamic coupling, as well as 
the interaction of the biopolymer with wall molecules 
in the region of the pore. Resolving these interactions 
is essential in understanding the physics underlying the 
translocation process. To this end, we model the dynam- 
ics of biopolymer translocation through narrow pores, us- 
ing the multiscale scheme described above. 

Our numerical simulations are performed in a three- 
dimensional box of size N x x N y x N z in units of the lattice 
spacing Ax. The box contains both the solvent and the 
polymer. We take N x = 2N y , N y = N z ; a separating 
wall is located in the mid-section of the x direction, at 
% = N x /2. For polymers with less than 500 beads we 
use N x — 80 while for larger polymers N x = 100. At 
t = the polymer resides entirely in the right chamber at 
x > N x /2. At the center of the separating wall, a square 
hole of side d — iAx is opened up, through which the 
polymer can translocate from one chamber to the other. 
A 3-D representation of a typical translocation event is 
shown in FigO The magnitude of the fluid speed is also 
mapped on different planes and as a 3-D contour surface 
surrounding the polymer beads. 

We elaborate next on the main parameters involved 
in the simulation (additional details are provided in Ref. 
Q). All parameters are measured in units of the lattice 
Boltzmann time-step At and spacing Ax, which are set 
equal to 1. The parameters for the Lennard- Jones poten- 
tial are a — 1.8, and e = 2 x 10~ 3 and the bond length 
among the beads is set to b — 1.2. The solvent density 
and kinematic viscosity are 1 and 0.1, respectively, and 
the inverse temperature is (3 = 10 . It must be noted 
that the friction coefficient taken as 7 = 0.1 is a parame- 
ter governing both the structural relation of the polymer 
towards equilibrium and the strength of the coupling with 
the surrounding fluid. The MD time-step is a 1/M frac- 
tion of the LB time-step, as mentioned previously, and 
we set M = 5. 

Translocation is induced by a constant electric force, 
Fdrive,i, which acts along the x direction and is confined 
in a rectangular channel of size 3Ax x Ax x Ax along the 
streamwise (x direction) and cross-flow (y,z directions). 
This force is included as an additional term in Eq. |T]) , and 
is the driving force representing the effect of the external 
field in the experiments. We use Fd r ive,i — 0.02. This 
choice of parameter values implies that we are describing 
the fast translocation regime, in which the translocation 
time is much smaller than the Zimm time, which is the 
typical relaxation time of the polymer towards its native 
(minimum energy, maximum entropy) configuration. Un- 
der these conditions, the many-body aspects of the poly- 
mer dynamics cannot be ignored because the beads along 
the chain do not move independently. 

Here, we model DNA as a polymeric chain of a number 
of segments (the beads) and trace its dynamic evolution 
interacting with a fluid solvent as it passes through a 
narrow hole that is comparable with the bead size. Each 
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FIG. 2: A typical translocation event: The 3-D box including the fluid and the polymer, as well as the wall (not shown) that 
separates the box into two chambers. A polymer (red beads) representing double-stranded DNA is shown at an instant where 
about 60% of the beads have translocated. The fluid is also shown as a 3-D representation of the velocity magnitude at the 
vicinity of the beads and on 2-D planes (red denotes high value). Translocation is induced by means of a pulling force Fdrive 
with a direction from right to left. 



bead maps to a number of base-pairs (bp) , ranging from 
about 8 (similar to the hydrated diameter of B-DNA in 
physiological conditions) to ~ 10 3 [2(|. In order to 
estimate this number for our simulations and interpret 
our results in terms of physical units we examine the 
persistence length (l p ) of the semiflexible polymers used 
in our simulations. We use the formula for the fixed- 
bond-angle model of a worm-like chain [2l| : 



lp 1 - cos(fl) 



(11) 



where (0) is complementary to the average bond angle 
between adjacent bonds. In lattice units (Ax) an average 
persistence length for the polymers considered was found 
to be approximately 12. For A-phage DNA, l p ~ 50 nm 
[22l |. which is set equal to l p for our polymers. Thereby, 
the lattice spacing is Ax ~ 4 nm, which is also the size 
of one bead. Given that the base-pair spacing is ~ 0.34 
nm, one bead maps to approximately 12 base pairs. With 
this mapping, the pore size is about ~ 12 nm, close to 
the experimental pores which are of the order of 10 nm. 
The polymers presented here with N = 20 — 500 beads 
correspond to DNA lengths in the range 0.2 — 6 x 10 3 bp. 
The DNA lengths used in the experiments are larger (up 
to ~ 10 5 bp); with appropriate computational resources, 
our multiscale scheme could handle these lengths. 

Having established the quantitative mapping of DNA 
base-pairs to the simulated beads, we seek a comparison 
of the statistical features of the simulated translocation 
process to experimental studies. The ensemble of our 
simulations is generated by different realizations of the 
initial polymer configuration to account for the statisti- 
cal nature of the process. We performed extensive sim- 
ulations of a large number of translocation events over 
100 — 1000 initial polymer configurations for each poly- 
mer length. The various events for a given length are de- 
picted in Fig. (3fa). The projected duration histograms 
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FIG. 3: (a) The variety of different translocation events and 
the corresponding translocation times; (b) from the events in 
(a) the duration histogram is extracted for a molecule with 
100 beads, which represents double-stranded DNA with 1.2 
kbp. Time is given in units of the LB time-step At. 



are shown in Fig. [3Jb) in LB units. Similar distributions 
were obtained for all the polymer lengths considered here, 
by accumulating all events for each length. By choos- 
ing lengths that match experimental data we compare 
the corresponding experimental duration histograms (see 
Fig. lc of Ref. [231) to the theoretical ones. This com- 
parison sets the LB time-step to At ~ 8 nsec. In Fig. Q] 
we show the time distributions for representative DNA 
lengths simulated here. In this figure, physical units are 
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used according to the mapping described above for di- 
rect comparison to similar experimental data [23]. The 
MD time-step for M = 5 will then be Atjvfo ~ 1.6 nsec 
indicating that the MD timescale related to the coarse- 
grained model that handles the DNA molecules is signif- 
icantly stretched over the physical process. Exact match 
to all the experimental parameters is of course not feasi- 
ble with coarse-grained simulations. Nevertheless, essen- 
tial features of DNA translocation are reasonably well 
reproduced, allowing the use of the current approach to 
model similar biophysical processes that involve biopoly- 
mers in solution. This can become more efficient by ex- 
ploiting the freedom of further fine-tuning the parameters 
used in the multiscalc model. 

The variety of all the different initial polymer realiza- 
tions produce a scaling law dependence of the transloca- 
tion times on length [l|| [23| . The duration histograms 
are not simple gaussians, but are rather skewed towards 
longer times. Accordingly, we use the most probable time 
(peak of the distribution shown by the arrow in Fig.[3Jb)) 
as the representative translocation time for every distri- 
bution; this is also the definition of the translocation time 
in experiments, to which we compare our results. Calcu- 
lating the most probable times for each length leads to 
the nonlinear relation between the translocation time r 
and the number of beads N: t(N) cx N a , with an ex- 
ponent a ~ 1.28 ± 0.01. The scaling law is reported in 
Fig. [5] and is in very good agreement with a recent ex- 
perimental study of double-stranded DNA translocation, 
that reported a ~ 1.27 ± 0.03 [23j. In the absence of 
a solvent the exponent rises to 1.36 ± 0.03. Such a dif- 
ference indicates a significant acceleration of the process 
due to hydrodynamic interactions. 

A. Dynamics of the translocation process 

We next turn to the dynamics of the biopolymer as it 
passes through the pore. The simulations confirm that 
the polymer moves through the pore in the form of two 
almost compact blobs on either side of the wall. One 
of the blobs (the untranslocated part, denoted by U) 
is contracting and the other (the translocated part, de- 
noted by T) is expanding. This behavior is visible in 
Fig. [5] for a random event and holds throughout the pro- 
cess, apart from the end points (initiation and completion 
of the translocation). A radius of gyration Ri(t) (with 
/ = U,T) can be assigned to each of these blobs, fol- 
lowing a static scaling law with the number of beads Ni: 
Riif) ~ Nj(t) with /i ~ 0.6 being the Flory exponent for 
a three-dimensional self-avoiding random walk. Based on 
the conservation of polymer length, Njj + Nt = N, an 
effective translocation radius can be defined as 

R E {t) = \R^> x {t)+R]i ti {t)f (12) 

which should be constant when the static scaling applies. 
We deduce from our simulations that -Re(£) is approxi- 
mately constant for all times throughout the process ex- 



cept near the end points, at which the polymer can no 
longer be represented as two uncorrelated compact blobs 
and the static scaling no longer holds. 

In Fig.[6][a), we represent the time evolution of all radii 
as averages over hundreds of events for a specific polymer 
length. The time shown is scaled so that t = 1 denotes 
the total translocation time for an event. By definition, 
Ru(t) vanishes at t — 1, while Rt increases monotoni- 
cally from t = up to t — 1, although it never reaches 
the value Rjj{t = 0) (we elaborate on this below). 

It is essential to check the validity of the static scaling 
with respect to the strength of the hydrodynamic field. 
To this end, the effective radii of gyration are further 
explored in relation to the parameter 7. We fixed the 
length to N = 200 beads and generated about 100 differ- 
ent initial configurations for each value of 7. We present 
the variation of Re with this coefficient in Fig. E^b): It 
is clearly visible in this figure that Re is almost constant 
for small 7 but as 7 increases, the radii are no more con- 
stant not only at the end points, but also throughout the 
translocation. Large values of the parameter 7 are inter- 
preted as a strong molecule-fluid coupling. The influence 
of the fluid on the beads, experienced by the back re- 
action, is large and suppresses the polymer fluctuations 
in such a way that the translocating biopolymer can no 
longer be represented as a pair of compact blobs, and the 
static scaling no longer holds. 

Inspection of all the biopolymers at the end of the event 
reveals that they become more compact after their pas- 
sage through the pore. This is quantitatively checked 
through the values of the radii of gyration: the radius of 
gyration is considerably smaller at the end than it was 
initially: R T (t = 1) < Ru(t = 0). The fact that as the 
polymer passes through the pore it becomes more com- 
pact than it was at the initial stage of the event may 
be related to incomplete relaxation, but this remains to 
be investigated. In Fig. [6]Jb), all effective radii of gy- 
ration at the final stage of the translocation decrease to 
a value smaller than the initial Rjj(t — 0). The ratio 
Rx{t — \)/Rjj(t — 0) is always smaller than 1 and ranges 
from 0.72 for 7=0. 1 to 0.90 for 7=0.5. 

Throughout its motion the polymer continuously in- 
teracts with the fluid environment. The forces that es- 
sentially control the process are the electric drive Fd r ive,i 
and the hydrodynamic drag Fd ra g,i which act on each 
bead. However, at the end points (initiation and com- 
pletion of the passage through the pore) entropic forces 
become important [24j | . The fluctuations experienced by 
the polymer due to the presence of the fluid are cor- 
related to these entropic forces which, at least close to 
equilibrium, can be expressed as the gradient of the free 
energy with respect to the fraction of translocated beads. 
At the final stage of a translocation event, the radius of 
the untranslocated part undergoes a visible deceleration 
(see Fig. (6fa)), the majority of the beads having already 
translocated. It is, thus, entropically more favorable to 
complete the passage through the hole rather than re- 
verting it, that is, the entropic forces cooperate with the 
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FIG. 4: Histograms of calculated translocation times for a large number of events and various DNA lengths shown in parentheses 
in kbp (10 3 bp). The arrows point to the most probable time for each polymer size. 
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FIG. 5: Scaling law of the translocation time with the DNA 
length. 



electric field and the translocation is accelerated. 

The entropic forces can also lead to rare events, such as 
retraction, which occur in our simulations at a rate less 
than 2% and depend on length, initial polymer configu- 
ration and parameter set. A retraction event is related 
to a polymer that anti-translocates after having partially 
passed through the pore. We have visually inspected the 
retraction events and associate them with the translo- 
cated part entering a low-entropy configuration (hairpin- 
like) subject to a strong entropic pull-back force from 
the untranslocated part: The translocated part of the 
polymer assumes an elongated conformation, which leads 
to an increase of the entropic force from the coiled, un- 
translocated part of the chain. As a result, the transloca- 
tion is delayed and eventually the polymer is retracted. 

The fact that the entropic forces are related to the 
number of translocated monomers Nx{t) led us to inves- 
tigate in more detail the time evolution of this quantity. 
The number of translocated monomers is plotted in Fig[7] 
for various initial configurations of the polymer. Each 
curve corresponds to a different completed translocation 
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FIG. 6: (a) Radii of gyration (translocated, untranslocated 
and effective) for a polymer representing 4.8 kbp as a function 
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are scaled with respect to the total translocation time and 
Ru(t = 0) for each case. 



event. The translocation for a given polymer proceeds 
along a curve closely related to its initial configuration 
and its interactions with the fluid; each polymer follows a 
distinct trajectory as indicated by the variety of curves. 
It is not possible to predict the ensuing behavior from 
the initial polymer configuration in a simple and unique 
manner. 
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FIG. 9: Distribution of the translocation work per unit time 
averaged over all events for a molecule with 3.6 kbp. The 
contributions from both the hydrodynamic and electric field 
are shown: Sfr(t) and respectively. 
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as functions of time for one translocation event of a polymer 
consisting of N = 300 beads (3.6 kbp). Time is scaled with 
respect to the total time of this event. 



B. Energetics of the translocation process 

As a final step, we study the work performed on 
the biopolymer throughout its translocation. On gen- 
eral grounds, hydrodynamic interactions are expected to 
minimize frictional effects and form a cooperative back- 
ground that assists the passage of the polymer through 
the pore. We investigate the cooperativity of the hydro- 
dynamic field through the synergy factor Sn(t) defined 
as the work made by the fluid on the polymer per unit 
time: 

Sli{t) = ^r = 7(f>w • *(*)) (is) 

i 

where Wjj is the work of the fluid on the polymer. 
Through this definition, positive values of this hydro- 
dynamic work rate indicate a cooperative effect of the 



solvent, while negative values indicate a competitive ef- 
fect by the solvent. The work done per timestep by the 
electric field (5b (t)) on the polymer can also be easily 
obtained through the expression: 

rlW N 
Mt) = ^=(Yl ' *<(*)) ( 14 ) 

i 

where We is the work of the electric drive on the poly- 
mer. The brackets in Eq. (fT3|) and (fT4)l denote averages 
over different realizations of the polymer for the same 
length. The results for the averages over all realizations 
are qualitatively similar to the work rates for an individ- 
ual event of the same length. For all lengths studied here, 
we found that the total work per timestep of the hydro- 
dynamic field (5ij(t)) on the whole chain is essentially 
constant, as shown in Fig. [5] for an individual event. For 
the same event, 5s(i) is also constant with time. In the 
same figure the kinetic energy K of the polymer (plot- 
ted as 7-K") for the same event is shown for comparison. 
The kinetic energy is also constant with time, as expected 
since the temperature in the simulations is held constant, 
but its fluctuations differ from those of Se and Sh- The 
larger value of K with respect to both Sh, Se can be 
justified by the fact that the bead velocities are larger 
compared to the fluid velocity. The hydrodynamic work 
per time is also larger than the corresponding electric 
field work, because the latter only acts in the small re- 
gion around the pore. The average of all these quantities 
over all events for the same polymer length are also con- 
stant, and show smaller fluctuations with time than those 
of any individual event. 

In addition to the variation of the work rates with time 
it is useful to analyze their distributions during transloca- 
tion events. We show these in Fig. \Q where it is evident 
that the distribution of Sh (t) lies entirely in the positive 
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range, indicating that hydrodynamics turns the solvent 
into a cooperative environment (it enhances the speed of 
the translocation process). In the same figure, the dis- 
tribution for Se (t) over all events for the same length is 
also shown; this distribution is mostly positive but has 
a small negative tail which indicates that beads can be 
found moving against the electric field. 

C. Performance data 

We turn next to some technical aspects of or multiscale 
simulations. The total cost of the computation scales 
roughly like 

t~(t LB V + t MD MN)N LB 

where £lb is the CPU time required to update a sin- 
gle LB site per timestep and Imd is the CPU time to 
update a single bead per timestep (including the over- 
head of LB-MD coupling); V is the volume of the com- 
putational domain in lattice units and N is the number 
of polymer beads, with M the LB-MD time-step ratio; 
Nlb is the number of LB timestcps. Due to the fact 
that the LB-MD coupling is local, the forces are short 
ranged and the SHAKE/RATTLE algorithms are empir- 
ically known to scale linearly with the number of con- 
straints, so that tMD is largely independent of N. The 
LB part is known to scale linearly with the volume oc- 
cupied by the solvent. Indeed, at constant volume, the 
CPU cost of the simulations scales linearly with the num- 
ber of beads: The execution times for 50, 100 and 400 
beads are 0.433, 0.489, and 0.882 sec/step, respectively 
on a 2GHz AMD Opteron processor. By excluding hy- 
drodynamics, these numbers become 0.039, 0.075, and 
0.318 sec/step. For the case where polymer concentra- 
tion is kept constant, the volume needed to accommo- 
date a polymer of N beads should scale approximately 
as A 18 . A typical translocation event with 500 beads, 
evolves over 30,000 LB steps or 150,000 MD steps. As- 
suming 250 flops/site/LB-step and 2500 flops/bead/MD- 
step, the previous equation leads to a computational cost 
of ~ 6 hrs. This is comparable with the time observed 
directly from the simulations (~ 7 hrs). 

IV. CONCLUSIONS 

We have presented a new multiscale methodology 
based on the direct coupling between atomistic motion 
and mesoscopic hydrodynamics of the surrounding sol- 
vent. Due to the particle- like nature of the mesoscopic 



lattice Boltzmann solver, this coupling proceeds via sim- 
ple interpolation/extrapolation in space and subcycling 
over time. Correlations between the atomistic and hy- 
drodynamic scales are also explicitly included through 
direct and local interactions between the solvent meso- 
molecules and the polymer molecules. As a result, hy- 
drodynamic interactions between the polymer and the 
surrounding fluid are explicitly taken into account, with 
no need of resorting to non-local representations, such as 
the long-range Oseen tensor used in Brownian dynam- 
ics. This allows a statc-of-thc-art modeling of biophysi- 
cal phenomena, where hydrodynamic correlations play a 
significant role. 

We have successfully applied our multiscale methodol- 
ogy to the problem of biopolymer translocation through 
nanoscale pores. Besides statistical properties, such as 
scaling exponents, the present methodology affords di- 
rect insights into the details of the dynamics as well as the 
energetics of the translocation process, thereby offering a 
very valuable complement to experimental investigations 
of these complex and fascinating biological phenomena. 
It also shows a significant potential to deal with problems 
that combine complex fluid motion and molecule dynam- 
ics. The efficiency of this scheme is also based on its rel- 
atively low computational demand. The molecular part 
is largely independent on the length of the molecule; the 
cost of the mesoscopic (LB) part is known to scale linearly 
with the volume occupied by the solvent. The linear scal- 
ing of the CPU time with the molecular size (at constant 
volume) is the key feature of the LB-MD approach, which 
permits the exploration of long biomolccules (N > 1000) 
with a relatively modest computational cost. Neverthe- 
less, resort to parallel computing is mandatory, and we 
expect the favorable properties of LB towards parallel 
implementations to greatly facilitate this task. This will 
also open the way to the simulation of systems at least an 
order of magnitude larger than those considered so far, 
making it possible to assign chemical specificity to the 
biopolymer constituents, rather than the generic nature 
of the beads that constitute the model polymers in the 
present study. Work along these lines is in progress. 
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